Emergent Criticality from Co-evolution in Random Boolean Networks 
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The co-evolution of network topology and dynamics is studied in an evolutionary Boolean network 
model that is a simple model of gene regulatory network. We find that a critical state emerges 
spontaneously resulting from interplay between topology and dynamics during the evolution. The 
final evolved state is shown to be independent of initial conditions. The network appears to be driven 
to a random Boolean network with uniform in-degree of two in the large network limit. However, 
for biologically realized network sizes, significant finite-size effects are observed including a broad 
in-degree distribution and an average in-degree connection between two and three. These results 
may be important for explaining properties of gene regulatory networks. 
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I. INTRODUCTION 

Boolean networks [1 B S H & H, 0, Hi have been 
extensively studied over the past three decades. They 
have applications as models of gene regulatory networks, 
and also as models of social and economic systems. As 
"coarse-grained" models of genetic networks they aim to 
capture much of the observed systematic behavior of the 
networks while simplifying local gene expression to a bi- 
nary (on/off) stateQ. Despite their simplicity, recent 
work has demonstrated that the model can indeed pre- 
dict the essential features of the dynamics of a biological 
genetic circuit [ToL ITlj|. An important feature of Boolean 
networks is that they have a continuous phase transi- 
tion between so-called ordered and chaotic phases. It 
has been argued [13, EH that gene regulatory networks 
of living systems should be at or close to criticality, at 
the so-called "edge of chaos" between the two phases, be- 
cause then they can maintain both the evolvability and 
stability. 

Many studies of critical random Boolean net- 
works (RBNs) have considered networks with homoge- 
neous topology in which each node has the same number 
of inputs from other nodes 0, EE EH E3, 
I2H |22| . Accumulating experimental evidence 
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however, shows that real genetic networks do not have ho- 
mogeneous connectivity, but, instead, are topologically 
heterogeneous. This diversity of architecture is, pre- 
sumably, of great importance for the stability of living 
cells. Studies of RBNs with heterogeneous topology have 
analytically determined the location of the ordered to 
chaotic phase transition in the large network limit and 
also demonstrated how different kinds of topology affects 
the stability of the dynamics [2rj I27I |28. 29). These stud- 
ies emphasized the importance of criticality and the in- 
fluence of the network's topology on its dynamics, but 
they did not attempt to explain how critical networks 
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with heterogeneous topology come to exist. 

Recent analysis of real gene regulatory networks (30l| 
has uncovered the possibility that the interactions be- 
tween genes can change in response to diverse stimuli. 
The resulting changes in the network topology can be far 
greater than what is expected simply from random mu- 
tation. However, few general principles are known about 
the evolution of network topology. In an effort to deter- 
mine what some of those principles may be Bornholdt 
and Rohlf[3l| studied a simple model of neural networks 
known as random threshold networks (RTNs). In their 
study, the topology of the RTN evolved according to a 
rule that depends on the local dynamics of the network. 
According to their rule, active nodes, whose binary state 
changes in time, tend to lose links, while inactive node, 
whose binary state is fixed, tend to gain new links. They 
observed that with this rule for changing network topol- 
ogy, in the limit of large networks, the RTNs evolved to 
a critical network with average connectivity K — 2. The 
study, therefore, discovered an interdependence between 
a network's dynamics and its topology. 

Motivated by these recent findings concerning gene 
regulatory networks and the behavior of RTNs, here we 
investigate the effect of a similar evolutionary rule on a 
simple model of a genetic regulatory system. In par- 
ticular, we study an evolutionary RBN model. How- 
ever, there is a fundamental difference between RTNs 
and RBNs. In RTNs the dynamics of each node is con- 
trolled by the same threshold function, while in RBNs the 
dynamics of each node is controlled by a randomly cho- 
sen Boolean function. Therefore, in addition to evolving 
the topology of the network, we also allow the Boolean 
functions used by the nodes to change. In the context of 
a model of gene regulation the rule we use to evolve the 
network topology assumes that there is some selection 
pressure on an individual gene due to its activity. The 
evolution causes genes that are in a frozen state of regu- 
lation, and, thus, are almost nonfunctional, to gain func- 
tionality, while it reduces the functionality of genes that 
are actively regulated. Thus this study investigates the 
co-evolution of network structure and network dynam- 
ics. We show that, independent of the initial topology of 
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FIG. 1: A directed graph G(7, 3) represents a homogeneous 
RBN with 7 nodes and in-degree connectivity of 3. Black and 
white represent binary states '1' and '0' respectively. The 
state vector of network is S(t) = (0, 1, 0, 0, 1, 1, 0). The arrow 
on each link indicates the direction of information flow in the 
network. 



setting the output value of the function for each set of 
input to be 1 with probability p. The bias p can be in- 
terpreted as a biochemical reaction parameter. 

Given the Boolean state of each node i at time t, 
(Ji(t), the state vector of network is defined as £(t) = 
(<7i(i), • • • , crjv(t)). The path that takes over time t 
is a dynamical trajectory in the phase space of system. 
Because the dynamics defined in Eq. [1] is deterministic 
and the phase space is finite, all dynamical trajectories 
eventually become periodic. That is, after some possi- 
ble transient behavior, each trajectory will repeat itself 
forming a cycle given by, 

£(t) = E(i + r). (2) 



the network, the RBN evolves to a critical network with 
a finite number of nodes and which has a heterogeneous 
topology. Two different variants of our model are consid- 
ered, and our principal conclusions are the same for both 
of them. Perhaps our most important result concerns 
the finite-size effects of the model. As the size of the 
network increases, the distribution of in-degree connec- 
tivity becomes increasingly narrow and sharply peaked at 
a value of K = 2. However, for biologically realized net- 
work sizes, we show that the final evolved critical state 
has a broad distribution of in-degree connectivity with 
an average value between 2 and 3. Both of these features 
also occur in real networks, suggesting that many of the 
topological features of real networks may be due to their 
finite size. 



The periodic part of the trajectory is the attractor of the 
dynamics, and the minimum V > that satisfies Eq. [2] is 
the period of the attractor. 

Two phases exist in RBNs, chaotic and ordered, char- 
acterized by their dynamical behavior 0, One impor- 
tant way of distinguishing the two phases is to measure 
the distribution of its attractor periods beginning with 
random initial states. For RBNs in the chaotic phase the 
distribution of attractor periods is sharply peaked near 
an average value that grows exponentially with system 
size N, and for RBNs in the ordered phase the distribu- 
tion of attractor periods is sharply peaked near an av- 
erage value that is nearly independent of N. Critical 
RBNs, however, have a broad power law distribution of 
attractor periods 14 1 . 



II. DEFINITION OF MODEL 

A. Dynamics of Random Boolean Networks 

A generalized RBN consists of N randomly intercon- 
nected nodes, i = 1,- ■■ ,N, each of which has Ki in- 
degree connections from nodes that regulate its behavior. 
The simplest Boolean network model is a homogeneous 
RBN in which each node has a same number of input 
nodes. In this case, the connections between nodes is de- 
scribed by a random directed graph G(N, K) consisting 
of N nodes with uniform in-degree K. Figure [T] illus- 
trates an example of G(7, 3). Each node has a Boolean 
dynamical state at time t, cr.;(i) = or 1. The state of 
each node at time t + 1 is a function of all states of its K , 
regulatory nodes at time t. Hence, the discrete dynamics 
of the network is given by 

<T i {t+l)=f i (<T il {t),a i2 (t),--- ,a iKi (t)) (1) 

where «i,i2,- ■ ■ ,£k 4 are those input nodes regulating node 
i. The function /, is a Boolean function of Ki variables 
that determines the output of node i for all of the 2 Ki 
possible sets of input. Note that random Boolean func- 
tions can be generated with an "interaction bias" p by 



B. Co-evolution in Random Boolean Networks 

In our model, the evolutionary changes of the topology 
of the network are driven by the dynamics of the network, 
and the functions that control the dynamics of network 
simultaneously evolve due to the changes in the network 
topology. Thus, this study investigates the co-evolution 
of network topology and dynamics. Similar to the one 
used in Ref. [3l| . the topology-evolving rule is simply 
that a frozen gene grows a link while an active gene loses 
a link. The dynamical functions can be changed in cither 
an annealed or a quenched way. The detailed algorithm 
is defined as follows: 

1. Start with a homogeneous RBN, G(N,Ko) with 
uniform in-degree connectivity Ki = Kq for all N, 
and generate a random Boolean function fa for each 
node i. 

2. Choose a random initial system state S(0). Up- 
date the state using Eq. Q] and find the dynamical 
attractor. See the appendix for a description of the 
algorithm used to find the attractor. 

3. Choose a node i at random and determine its av- 



erage activity 0(i) over the attractor. 

T+r-i 

0(i) = ± E ( 3 ) 

where T is a time large enough so that the periodic 
attractor has been reached, and T is the period of 
the attractor. If 0(i) = 1 or 0, then its state does 
not change over the duration of the attractor; it is 
frozen. Alternatively, if < 0(i) < 1, then node i 
is active during the attractor. 

4. Change the network topology by rewiring the con- 
nections to the node chosen in the previous step. If 
it is frozen, then a new incoming link from a ran- 
domly selected node j is added to it. If it is active, 
then one of its existing links is randomly selected 
and removed. Note that this rewiring changes Ki. 

5. The Boolean functions of network are regenerated. 
Two different methods have been used: 

• Annealed model: A new Boolean function is 
generated for every node of the network. 

• Quenched model: A new Boolean function is 
generated only for the chosen node i, while the 
others remain what they were previously. 

6. Return to step 2. 

The time scale for an evolutionary change of the net- 
works, steps 2 — 6 above, is called an epoch. As we will 
see, using this rewiring rule the network topology evolves 
from a homogeneous one to a heterogeneous one. For 
simplicity, all random Boolean functions are generated 
with p = 1/2, and therefore all Boolean functions with 
the same in-degree are equally likely to be generated. 

III. SIMULATION AND RESULTS 

We have simulated both the annealed and the 
quenched variants of the model. Both variants give very 
similar results, and our principal findings are the same 
for both variants. Therefore, we will present here mainly 
results from the annealed variant. Graph (a) of Fig. [3] 
shows the evolution of the average in-degree connectiv- 
ity K = -i? Ki for networks of size N = 30 in the 
annealed variant of the model. Four curves are shown. 
They show the results obtained by beginning with net- 
works with different uniform connectivity Kq = 2, 3, 4, 
and 5. Each curve is the average of 15,000 independent 
realizations of the network evolution. This ensemble av- 
erage is indicated by the angular brackets. Each different 
realization in an ensemble begins with a different random 
network and with a different random initial state vec- 
tor. Remarkably, despite the difference in initial condi- 
tions, all four curves collapse after about 200 epochs, and 
they all approach the same final statistical steady state 
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FIG. 2: (color online) (a). Evolution of the ensemble averaged 
in-degree connectivity in the annealed model for networks of 
size N = 30. The networks in each ensemble initially start 
from different uniform connectivity, Ko — 2, 3, 4, and 5, but 
reach a same statistical steady state (K) = 3.06. Each ensem- 
ble contained 15,000 realizations of the network, (b). Evolu- 
tion of ensemble averaged in-degree connectivity for networks 
of three different size N — 30, 50, and 100 in the annealed 
model. 



that has an average in-degree connectivity {K) — 3.06. 
This occurs without tuning and suggests that the final 
evolved topology of the network is independent of the 
initial topology of the network. The steady state value 
of (K) depends on the size of the system as shown in 
graph (b) of Fig. [21 Starting with networks that all have 
the same initial uniform connectivity Kq = 4, but which 
have different size N = 30, 50, and 100, we find that 
larger networks evolve to steady states with smaller val- 
ues of (K). Very similar results are obtained for the 
quenched version of the model. For example, for net- 
works with N = 30 the steady state value of the average 
connectivity is (K) = 3.08. 

We have also calculated the in-degree and out-degree 
connectivity distributions, P{Ki n ) and P(K out ), of the 
evolved RBNs in the steady state. Initially, all nodes 
of the network have a uniform in-degree Ko, meaning 
that the in-degree distribution is a discrete delta func- 
tion P(Ki n ) — Sx in ,K , and the out-degree distribution 
is a binomial distribution. However, through the evolu- 
tionary rewiring of the network both the in-degree and 
out-degree distributions change. The in-degree and out- 
degree distributions in the steady state of the annealed 
version, are shown in Fig. [3] for N — 200. They are both 
right skewed bell-shaped distributions peaking at K = 2. 
The out-degree distribution remains a binomial distribu- 
tion but the average connectivity changes. The in-degree 
distribution, although it has the same average connectiv- 
ity as the out-degree distribution, is more sharply peaked. 
As the size of the network grows, the in-degree distribu- 
tion becomes increasingly narrow and peaked at the value 
Ki n — 2, as shown in Fig. 2J Based on this observation, 



5 10 15 

K 

FIG. 3: (color online) Distribution of in-degree (square) and 
out-degree (circle) connectivity in the annealed model. The 
size of the networks is N = 200. 
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FIG. 5: (color online) Power law distribution of steady state 
attractor period T in both annealed (circle) and quenched 
(square) models for N = 200 system. The dashed straight 
line has a slope of 1.0. 




FIG. 4: (color online) Distribution of in-degree connectivity 
in the annealed model. The network sizes are N = 200, 400, 
and 1000. 



evolutionary RBNs. 

Given the steady state value (K) = 2 in the large net- 
work limit N — > oo, we studied the finite-size effects in 
the model. As shown in Fig.[6j the values of (K(N)) for 
finite N obey the scaling function 

(K(N)) -2 = AN' 13 . (4) 

Fitting the data to this function, we find that the co- 
efficient is A = 2.50 ± 0.06 and the exponent is [5 = 
0.264 ±0.005. Thus the value of (K(N)) is always larger 
than 2 for finite N. Note that steady state values of the 
average connectivity in random threshold networks have 
a similar scaling form [3l|, but, in that case A — 12.4±0.5 
and = 0.47 ±0.01. 

IV. DISCUSSION AND CONCLUSIONS 



we conjecture that the distribution tends to converge into 
a discrete delta function <5jQ n ^ in the large network limit 
N — > oo, indicating that the network becomes a homo- 
geneous RBN in that limit. 

In order to probe the dynamical nature of evolved 
steady states we computed the distribution P(T) of 
steady state attractor period T in the ensemble of RBNs 
simulated. The distribution has a broad, power-law be- 
havior for both the annealed and quenched variants of 
the model. Figure shows the results for networks with 
N = 200. As long as N is about 30 or larger, results 
for other size networks are similar. As discussed above, 
this power-law distribution indicates that the networks 
have critical dynamics. Also in the figure, the straight 
line has a slope of 1.0. Thus, the critical exponent de- 
scribing the power-law is approximately 1.0. This value 
of the exponent is obtained for all system sizes studied in 
both the quenched and annealed versions of the model. 
In short, we find that a robust criticality emerges in the 



The mechanism we find here that leads to the emer- 
gent critical state has some similarity to self-organized 
criticality (SOC)0, but is different. SOC is the 
tendency of driven dissipative dynamical systems to orga- 
nize themselves into a critical state far from equilibrium 
through avalanches of activity of all sizes. In our particu- 
lar model, the evolutionary Boolean network is a dissipa- 
tive dynamical system because multiple different states 
may map into the same attractor so that information is 
lost. Similar to SOC systems, our model is driven subject 
to two competing rules, and the network organizes itself 
into a steady state that results from a dynamical balance 
of the competition between those rules. Moreover, the 
critical state is robust irrespective of initial connectivity 
in both the quenched or annealed versions of the model. 
The emergent critical state acts like a global attractor in 
the evolution process. However, unlike the mechanism 
of traditional SOC, but similar to the mechanisms that 
have been shown to lead to criticality in random thresh- 
old networks (3lT | and in homogeneous Boolean networks 
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FIG. 6: (color online) Finite-size effects in the annealed 
model. The data shown are for systems of different sizes 
N = 30, 50, 80, 100, 150, 200, and 400. The dashed straight 
line has a slop of 0.26. 



[201 ] . the self-organizing mechanism here is based on a 
topological phase transition in dynamical networks. 

Our results indicate that, with the rewiring rules we 
use, networks in the limit N — > oo will evolve to have a 
homogeneous in-degree connectivity of 2. However, we 
find that finite size networks evolve to have a broadly 
distributed heterogeneous in-degree connectivity. The 
average in-degree connectivity of the evolved networks 
is between 2 and 3 for biologically realized network sizes. 
This result may be important for explaining the observed 
structure of real gene regulatory networks. Real genetic 
networks have a number of nodes N ranging from near 
100 to thousands, and typically have a heterogeneous 
connectivity with an average in-degree slightly larger 
than 2. For example, a recent experiment studying the 
gene regulatory network of S. cerevisiae [30( found that 
the network contains 3420 genes and has an average in- 
degree connectivity Ki n =2.1. 

Finally, we note that real genetic networks exhibit 
an approximately scale-free out-degree distribution while 
the in-degree distribution is exponentially decaying[30, 
34]. In our numerical study, we similarly obtain an in- 
degree distribution that decays faster than the out-degree 
distribution, but our model does not produce a scale-free 
like out-degree distribution. Therefore, in order to be 
more realistic the model needs to be extended by adding 
other factors that will capture this feature. We do not 
believe though that such extensions will alter the princi- 
pal conclusions of this paper concerning the importance 
of finite size effects in the evolution of network topology 



in real genetic networks. 
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APPENDIX 

In the work presented here and in our previous 
studies[2(| [U [22j the following algorithm was used to 
determine the dynamical attractor. 

Begin by using the initial state 53(0) as the "check- 
point" state. Then for each time T, < T < Ti, the 
state is updated using Eq. 1. If 53 (T) = 53(0), then the 
attractor is found, it has period T, and the search ends. 
Note that the new state is compared only to the check- 
point state and no other previous states. 

If after T\ updates no attractor is found then the 
checkpoint state is changed to 53 (Ti). For each time T, 
Ti < T < T2, the state is again updated using Eq. 1. If 
53(T) = 53(Ti), then the attractor is found, it has period 
T — Ti, and the search ends. 

If after T 2 updates the attractor has still not been 
found then the checkpoint state is changed to 53 (Ti). For 
each time T, T 2 < T < Ti + T max , the state is again up- 
dated using Eq. 1. If 53 (T) = 53 (T 2 ), then the attractor 
is found, it has period T — T 2 , and the search ends. 

Finally, if after T 2 + T max updates the attractor has 
still not been found then the state 53 {Ti + T max ) is used 
as the new checkpoint state. If 53 (T) = 53 (T 2 + T max ), 
then the attractor is found, it has period T — Ti — T max , 
and the search ends. 

If no attractor is found after this procedure, then we 
stop. In this case, the average output state in Eq. 3 
is calculated assuming that the attractor period is T — 
Tmax and that the final checkpoint state 53 (Ta + T max ) 
is on the attractor. 

This algorithm finds all attractors that have period less 
than or equal to T max and that have a transient time 
to reach the attractor from the initial state less than or 
equal to Ti + T max . For the results presented in this 
paper we used Ti = 100, T 2 = 1000, and T max = 100000. 



[1] S. A. Kauffman, J. Theor. Biol. 22, 437 (1969). 

[2] S. A. Kauffman, Physica D 42, 135 (1990). 

[3] B. Derrida and Y. Pomeau, EuroPhys. Lett. 1, 45 (1986) 

[4] B. Derrida and G. Weisbuch, J. Phys. 47, 1297 (1986). 

[5] U. Bastolla and G. Parisi, Physica D 115, 203 (1998). 



[6] R. Albert and A. L. Barabasi, Phys. Rev. Lett. 84, 5660 
(2000). 

[7] M. Aldana, S. Coppersmith, and L. P. Kadanoff, in Per- 
spectives and Problems in Nonlinear Science, edited by E. 
Kaplan, JE. Marsden, and K.R. Sreenivasan (Springer- 



6 



Verleg, Berlin, 2003). 
[8] B. Luque, F. J. Ballesteros, and E. M. Muro, Phys. Rev. 

E 63, 051913 (2001). 
[9] S. Bornholdt, Science 310, 449 (2005). 
[10] R. Albert and H. G. Othmer, J. Theor. Biol. 223, 1 

(2003). 

[11] F. Li, T. Long, Y. Lu, Q. Quyang, and C. Tang, Proc. 
Natl. Acad. Sci. U.S.A. 101, 4781 (2004). 

[12] N. H. Packard, in Dynamics Patterns in Complex Sys- 
tems, edited by J. A. S. Kelso, A. J. Mandell, and M. F. 
Shlesinger( World Scientific, Singapore, 1988). 

[13] S. A. Kauffman, The Origins of Order (Oxford Univer- 
sity Press, New York, 1993). 

[14] U. Bastolla and G. Parisi, J. Theor. Biol 187, 117 (1997). 

[15] F. Greil and B. Drossel, Phys. Rev. Lett. 95, 048701 
(2005). 

[16] B. Drossel, Phys. Rev. E 72, 016110 (2005). 

[17] B. Drossel, T. Mihaljev, and F. Greil, Phys. Rev. Lett. 

94, 088701 (2005). 
[18] B. Samuelsson and C. Troein, Phys. Rev. Lett. 90, 

098701 (2003). 

[19] J. E. S. Socolar and S. A. Kauffman, Phys. Rev. Lett. 

90, 068702 (2003). 
[20] M. Paczuski, K. E. Bassler, and A. Corral, Phys. Rev. 

Lett. 84, 3185 (2000). 
[21] K. E. Bassler, C. Lee, and Y. Lee, Phys. Rev. Lett. 93, 



038101 (2004). 

[22] K. E. Bassler and M. Liu, in Noise in Complex Systems 

and Stochastic Dynamics III, Proc. of SPIE Vol. 5845, 

edited by L. B. Kish, K. Lindenberg, and Z. Gingl (SPIE, 

Bellingham, WA, 2005), P. 104. 
[23] D. Thieffry, A. M. Huerta, E. Perez-Rueda, and 

J. Collado-Vides, Bioessays 20, 433 (1998). 
[24] T. I. Lee, N. J. Rinaldi, F. Robert, D. T. Odom, and 

Z. B. J. et al, Science 298, 799 (2002). 
[25] A. H. Y. Tong, G. Lesage, G. D. Bader, H. Ding, and 

H. X. et al, Science 303, 808 (2004). 
[26] B. Luque and R. V. Sole, Phys. Rev. E 55, 257 (1997). 
[27] J. J. Fox and C. C. Hill, Chaos 11, 809 (2001). 
[28] M. Aldana, Physica D 185, 45 (2003). 
[29] M. Aldana and P. Cluzel, Proc. Natl. Acad. Sci. U.S.A 

100, 8710 (2003). 
[30] N. M. Luscombe, M. M. Babu, H. Yu, M. Snyder, and 

S. A. T. M. Gerstein, Nature 431, 308 (2004). 
[31] S. Bornholdt and T. Rohlf, Phys. Rev. Lett. 84, 6114 

(2000). 

[32] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 
59, 381 (1987). 

[33] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. A 38, 
364 (1988). 

[34] R. Albert, J. Cell Sci. 118, 4947 (2005). 



